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Abstract. We have used the master equation approach to study a moderately complex network of diffusive 
reactions occurring on the surfaces of interstellar dust particles. This network is meant to apply to dense clouds 
in which a large portion of the gas-phase carbon has already been converted to carbon monoxide. Hydrogen 
atoms, oxygen atoms, and CO molecules are allowed to accrete onto dust particles and their chemistry is followed. 
The stable molecules produced are oxygen, hydrogen, water, carbon dioxide (CO2), formaldehyde (H2CO), and 
methanol (CH3OH). The surface abundances calculated via the master equation approach are in good agreement 
with those obtained via a Monte Carlo method but can differ considerably from those obtained with standard 
rate equations. 
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1. Introduction 

Rate equations have been widely used in simulations of 
gas-phase processes in the interstellar medium (Le Teuff 
et al. 120000 . Their application has also been extended to 
treat diffusive reactions on the surfaces of dust particles 
(Pickles & Williams IHT77I Hasegawa et al. HM2|) . This ap- 
proach, however, is valid only if the average number of 
reactive surface species per grain is large and the discrete 
nature of the system can be neglected. When this num- 
ber becomes small, the rate equations may no longer con- 
stitute an accurate description of the chemical processes 
(Tielens & Haeen IT982fl : moreover, for species with sur- 
face populations less than one, it can become meaningless 
to use the rate equations (Charnley et al. 119971 Caselli et 
al. lTMgjl . 

This problem has spurred attempts to develop alterna- 
tive methods for handling diffusive reactions on granular 
surfaces. One possible solution is to use Monte Carlo pro- 
cedures to simulate the grain-surface chemistry. Different 
methods, based on Monte Carlo simulations, have already 
been employed for various grain surface networks (Tielens 
& Hagen 113521 Charnley et al. ITW71 Charnley I2UUT1 
Caselli et al. 2002). These simulations were performed un- 
der the constraint that, during the evolution of the surface 
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chemistry, the gas-phase abundances of all species stay 
constant - a requirement that cannot be met in complex 
gas-grain models, where it is essential that the gas-phase 
and the grain-surface reactions run in a parallel mode dur- 
ing the calculations and that both adsorption onto and 
desorption from grain surfaces occur (Ruffle & Herbst 
2000). ft does not appear possible, however, to follow si- 
multaneously the gas-phase chemistry with rate equations 
and the surface chemistry with a Monte Carlo approach. 
Although one can use the Monte Carlo method to simulate 
gas-phase reactions (Charnley 1998), a Monte Carlo solu- 
tion for both gas-phase and grain-surface chemisty could 
be done only with an exceedingly large amount of comput- 
ing time, and has not yet been attempted. Current gas- 
grain chemical models (e.g. Ruffle & Herbst 20QUJ use a 
semi-empirical modification of the standard rate equations 
for diffusive surface chemistry known as the " modified rate 
approach" (Caselli et al. |P9"9"%1 Stantcheva et al. I2OTT1 
Caselli et al. 2002). Although this approach is efficient, its 
semi-empirical nature raises doubts of its suitability under 
all conditions. 

Recently, two groups (Biham et al. 120011 Green et 
al. 12)01) proposed a second stochastic approach to granu- 
lar chemistry, known as the master equation treatment. In 
this approach, differential rate equations for species with a 
small surface abundance are replaced by differential equa- 
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number of atoms or molecules of that species (0,1,. .. ) is 
present on a grain at any time. In general, the probabili- 
ties for each surface species are not independent and so an 
exact treatment requires the determination of joint proba- 
bilities (e.g., for of species A, 1 of species B, etc.). Green 
et al. I|2001[l used the master equation method to study 
the simple O, H system; in this system O and H atoms, 
with fixed gas-phase abundances, are allowed to adsorb 
onto grains and react to produce the three molecules OH, 
O2, and H 2 via diffusive processes. For this system, de- 
termination of the surface populations of the minor but 
important species H and O requires the calculation of 
two-body probabilities, while surface abundances of the 
diatomic species can be calculated from normal rate equa- 
tions. The two-body probability approach to the simple O, 
H system can itself be simplified by a number of approx- 
imate methods (Green et al. |2001). One approximation, 
proposed by Biham et al. H2001fl but not attempted by 
Green et al. Ij2001(l . reduces the two-body probability to 
independent probabilities for the individual species O and 
H. With this approach, the so-called "many-body" mas- 
ter equation reduces to separate master equations for the 
individual species. Although the approximation leads to 
a significantly smaller number of simultaneous differential 
equations to solve, its validity is not clear for systems with 
strong correlations between the surface abundances of the 
different minor species. 

If one attempts to scale up the many-body master 
equation approach to model realistic complex networks of 
surface reactions, one can experience serious problems in- 
volving both computer time and memory, unless suitable 
approximations can be found. In order to learn more about 
the use of the method for larger systems, we have con- 
sidered an intermediate system more complex than previ- 
ously studied but not at the level of complexity needed for 
a complete network of surface reactions. In particular, we 
report here the use of the many-body master equation ap- 
proach to solve a grain-surface chemical network in which 
gaseous H, O, and CO accrete onto grain surfaces and lead 
to the production of molecular hydrogen, formaldehyde, 
methanol, water, oxygen, and carbon dioxide (Caselli et 
al. 120021 see also Charnley et al. IT9T)7|. This system has 
previously been studied by the modified rate and Monte 
Carlo approaches (Caselli et al. I2002[) . We consider here 
a variety of different diffusive rates and temperatures to 
see how the master equation approach fares. As opposed 
to the previous treatments based on the master equation, 
we utilize a time-dependent approach designed eventually 
to be coupled with a time-dependent gas-phase chemistry, 
since the advantage of the direct solution of the master 
equation compared with the Monte Carlo realization is 
that the differential equations are easily coupled to rate 
equations for the gas-phase species. 

The organization of the paper is as follows. In the next 
section, we discuss the chemical network of diffusive sur- 
face reactions and the different rates chosen. In Sect. 3, we 
write out the differential equations needed to model the 



Table 1. H, O, and CO gas-phase abundances (cm 3 ) 
utilized 



Abundance n 


Low 


Intermediate 


High 


H 


1.15 


1.15 


1.10 





0.09 


0.75 


7.0 


CO 


0.075 


0.75 


7.5 



in Sect. 4 we present our results, and compare them with 
Monte Carlo and rate approaches. A general discussion is 
contained in Sect. 5. In addition, a detailed mathematical 
discussion of the master equation and Monte Carlo realiza- 
tions to stochastic kinetics is contained in the Appendix. 

2. The H, O, CO Network 

In this network, three gas-phase species - H, O, and CO 
- are allowed to accrete onto a granular surface where 
they react via diffusion to produce the stable molecules 
H 2 , 2 , H 2 0, H 2 CO, CH3OH, and C0 2 , as well as reac- 
tive intermediate species. Of the ten reactions, nine are 
association reactions, in which a single product is formed. 
No gas-phase chemistry is allowed to occur and the gas- 
phase concentrations of the three species are assumed to 
remain constant despite accretion onto grains. This some- 
what contradictory assumption permits the surface chem- 
istry to occur only over a period sufficiently short that 
the abundances of the gas-phase species do not change 
dramatically. 

The calculations were carried out for three different 
sets of gas-phase abundances of the accreting species, 
known respectively as the low, intermediate, and high den- 
sity cases. All of these cases refer to dense clouds in which 
most of the atomic hydrogen has already been converted 
into molecular hydrogen, leaving only a small remnant in 
atomic form. Such conditions pertain when a significant 
amount of CO has been produced via gas-phase chem- 
istry. The abundances n(cm -3 ) of H, O, and CO shown in 
Tabled were obtained from steady-state gas-phase mod- 
els run at molecular hydrogen densities of 10 3 , 10 4 , and 
10 5 cm~ 3 . In the low-density case, there is more atomic 
hydrogen than O and CO, and the chemistry will be seen 
to be strongly reductive, whereas in the high-density case, 
there is little atomic hydrogen around to react with CO. 

The surface reactions for the network are listed in 
Table along with activation energies in K where ap- 
propriate. The activation energies E a (K) are approxi- 
mate only (Caselli et al. I2002j) . The total number of sur- 
face species is 12; this includes the highly reactive radi- 
cals OH, HCO, and H3CO. The key reaction sequence in 
the network is the relatively slow hydrogenation of CO 
into methanol (CH3OH) via H-atom addition reactions. 
Unambiguous laboratory evidence for this hydrogenation 
is not available (Hiraoka et al. 120001) . but conditions in 
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Table 2. Surface reactions in the H,0,CO model 



Number Reaction E a (K) 



1 


H 


+ H 2 




2 


H 


+ O — ► OH 




3 


H 


+ OH — ► H 2 




4 


H 


+ CO — ► HCO 


2000 


5 


H 


+ HCO — ► H 2 CO 




6 


H 


+ H2CO — > H3CO 


2000 


7 


H 


+ H3CO — ► CH3OH 




8 





+ — > 2 




9 





+ CO — > co 2 


1000 


10 





+ HCO — ► CO2+H 





Table 3. Energy barriers against diffusion (low values) 
and desorption energies (K) 



Species 


E h (K) 


E v (K) 


H 


100 


350 


O 


240 


800 


OH 


378 


1260 


H 2 


135 


450 


2 


363 


1210 


H2O 


558 


1860 


CO 


363 


1210 


HCO 


453 


1510 


H2CO 


528 


1760 


CH3O 


651 


2170 


CH3OH 


618 


2060 


CO2 


750 


2500 



complete synthesis in Tabled is by no means ruled out by 
experiments. Once produced, all stable species except H 2 
remain on the grain surface; the evaporation of molecular 
hydrogen is included. 

Whether one uses rate equations, the Monte Carlo ap- 
proach, or the direct master equation method, it is neces- 
sary to utilize diffusion rate coefficients k for the reactive 
surface species (Hasegawa et al. 1992 Appendix). The 
rate coefficients here are in units of s , as preferred by 
Casern, Hasegawa, & Herbst (1998). These are the sum of 
the rates (s _1 ) of the reactive partners to traverse an 
entire grain, which is here assumed to contain 10 6 binding 
sites, multiplied by a factor k that accounts for any non- 
zero chemical activation barrier (Hasegawa et al. I1992[l . 
The rates depend strongly on the barriers against diffu- 
sion Eh from site to site chosen, and whether diffusion 
occurs via thermal hopping or via quantum mechanical 
tunneling (Tielens & Haeen 11352)) . 

In our calculation, we have considered three sets of 
barriers against diffusion. The first, which comes from 
the earlier astrochemical literature (Allen & Robinson 
113771 Tielens & Hagen[]JH2 Hasegawa et al. con- 
tains rather low barriers and allows efficient tunneling for 
atomic H. The second and third are based on the recent 
experiments of Pirronello et al. I|19971 Q')99) as simulated 
by Katz et al. 1|1999[) . which show that atomic H moves 
much more slowly on olivine and amorphous carbon than 
previously assumed in the astrochemical literature. Two 
sets of barriers based on these experiments on olivine have 
been used (Ruffle & Herbst 200U) - one, designated "slow 
H", in which only the H atom barrier is raised, and the 
other, designated "slow", in which all other barriers are 
raised proportionately. In both sets of barriers, no tun- 
neling is allowed, since no tunneling of H was detected in 
the laboratory. For our calculations here, the earlier as- 
trochemical values are used principally because the slower 
diffusion rates cannot produce much formaldehyde and 
methanol in the small times considered, and so are not 
emphasized (Caselli et al. 2002 ). In our more complex gas- 



Table 4. Assorted rates for selected species at 10 K 



Species 


fcacc (cm 3 S 1 ) 


^evap (s ) 


(s" 1 ) 


H 


1.45(-5) 


1.88(-3) 


5.14(+4) a 


O 


3.62(-6) 


2.03(-23) 


4.24(-5) 


CO 


2.73(-6) 







a Quantum tunelling included 



that formaldehyde and methanol can be produced over 
long periods of time (Ruffle & Herbst [2000 ) . 

In addition to diffusive rates, the rates of adsorp- 
tion and desorption must be included in our calculation 
(Hasegawa et al. H992fl . Adsorption is assumed to occur at 
unit efficiency once a gas-phase species strikes a grain. We 
consider only thermal desorption (evaporation) and treat 
it as in previous models (see e.g. Caselli et al. 199.8); the 
rates, exponentially dependent on the desorption energy 
-Ed) are included for the accreting gas-phase species and 
for molecular hydrogen product. Heavy molecular species 
desorb too slowly for this process to be considered here. 
The small barriers against diffusion and the desorption 
energies for all species in the model are listed in Table |3J 

Table 0] gives our values for the accretion rate coeffi- 
cients fc acc (cm 3 s _1 ) onto a grain, the evaporation rates 
^cvap from the grain, and the diffusion rates for the species 
H, 0, and CO at 10 K, unless they are vanishingly small. 
Calculations have been done mainly at this temperature, 
although temperatures up to 20 K have been considered. 

3. Master equation for H,0,CO system 

Of the three species - H, O, and CO - that accrete onto 
grain surfaces, the first two are very reactive and never 
build up large surface populations. On the other hand, 
CO reacts only slowly via reactions with activation en- 
ergy and so can build up a large surface population under 
certain circumstances. In general, all surface species in our 
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The major ones correspond to the stable species that re- 
act slowly if at all (CO, H 2 CO, CH 3 OH, H 2 0, 2 , H 2 , 
C0 2 ) and can build up large abundances, while the mi- 
nor ones are atoms and radicals likely to have a surface 
number, defined as the number of species per grain, at or 
below unity. We treat the minor species - H, O, OH, HCO, 
and H3CO - probabilistically with corresponding surface 
numbers i\, £2, £3, £4, and £5. 

The first step in the master equation approach is to 
solve for the joint probability -P(ji,«2)*3,*4,*5), defined as 
the probability that i\, £2, £3, £4, and £5 numbers of minor 
species exist on the surface as a function of time. In the 
calculations discussed here, we start with the initial con- 
dition that the joint probability is unity for P(0, 0, 0, 0, 0). 
Let X and Y be, respectively, the ] th and k th reactive ele- 
ments of the ordered set {H, O, OH, HCO, H 3 CO }. Let 
Z represent any of the major species. The time derivative 
of the five-body probability for each value of £ 1 , £ 2 , £3 , £4 , £5 
can then be written as (Appendix) 



5]fc acc (X)n(X) [P(. 

{X} 



1,...)-P(.. 



•)] 



-^-i^X) [fe + l)P(...,ij + 1, ...) - ';/'!..../,. ...) 

{X} 

- £ fc x , Y (£j + l)(£)fc + l)P(...,«j + l ) - ! H + l,-) 

{X,Y} 

- ^2 kx,Y(ij)(ik)P(-,ij, -,h,-) 

{X,Y} 



(1) 



{X} 
{X} 



(i j +2)(i j + l) 



P(...,£ J +2,...) 



P(..,£i,...) 



]T (JVz>*x,z[(i 3 - + l)P(...,^ + l,...)] 



{X,Z} 



E (^z) fc x,z [ijP(. 



1, j, . 



{X,Z} 



where n stands for gas-phase concentration and (N) for 
surface abundance; i.e., the average number of atoms or 
molecules per grain of a species. 

The first term on the right-hand side of eq. Q ac- 
counts for changes in the state of the surface on a partic- 
ular grain due to accretion of species. In this particular 
case, the sum consists of two terms because the only ac- 
creting species with minor surface abundance are H and 
O. The second term describes the changes of the system 
due to evaporation, and the remaining terms take into ac- 



are subdivided into expressions for reaction between two 
distinct minor species (X,Y), for self-reaction (X,X), and 
for reaction between a minor and a major species (X,Z). 
Note that all of the terms refer to reactants; there is also 
one minor product - H atoms in reaction 10 of Table |2J 
To include the production of H in Eq. Q requires a term 
which, when X and Y are equal to O and HCO, contains 
probability functions where three indices change. We have 
not included this term in eq. Q for simplicity, but it is 
of course included in our calculations. The average abun- 
dances of the major species are obtained from rate equa- 
tions discussed below. It is easily shown that the total 
probability as a function of time is conserved at unity. 

Because the abundances of minor species are low, joint 
probabilities with high numbers of these particles are very 
unlikely and therefore, the probabilities for such states 
can be neglected. In particular, we choose a set of pa- 
rameters Af = {Afi, A/2, A3, A/4, A/5 } such that the only 
probabilities P(...,ij,...) to be considered possess ij < A/}, 
je {1,2,3,4,5}. Specific choices for the set Af are dis- 
cussed in the Results section. It is obvious that for the 
master equation method to be feasible, the set Af must 
contain elements as small as possible. 

Once the probabilities are determined by integration 
for a specific time, the average numbers of minor species, 
(-ZVh), (No), etc., as well as the correlation terms (NuNo), 
(NhNqh), etc., can be calculated from these probabilities; 



e-g-, 



hP(h,h,h,ii, £5)- 



(2) 



Both average numbers and correlations are then used in 
the rate equations for the abundance of the major surface 
species, while the latter can also be used to test how in- 
dependent or correlated the minor species are. In the rate 
equations for major species, shown below, the division of 
species into major and minor ones leads to the fact that 
correlations are only used for pairs of minor species: 



=-t-i p (H 2 ) (Y H2 ) 
+ fc H ,H x 0.5 x (N H {N H - 1)>, 



d(Np 2 ) 
dt 

dt 



= k Q ,o x 0.5 x (N (N - 1)}, 



"H.OH 



X (N h N h), 



d(Nco) _ 



~3t ~ ' 



3 (CO)n(CO) 

d(N K2 co) 



ko,co(N C o){N ) 



"H.HCO 



{NnN: 



hco; 



(3) 

(4) 
(5) 

(6) 
(7) 
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•5 



d{N { 



CH 3 OH/ 



dt 



= fcH,H 3 CO (NaNn 3 co) 



(8) 



d{N £° 2) = ko,co(N C o)(No) + k OM co{N N W o)- ( 9 ) 



3.1. Approximation of Independent Probabilities 

Considering the large number of coupled differential equa- 
tions for the many-body P required in the exact master 
equation approach, it is perhaps useful to consider whether 
the use of one-particle probability functions is adequate. 
This approach was suggested but not tested by Biham 
et al. H2001f) f° r the simple O, H system, and leads to 
a different and somewhat simpler system of differential 
equations. The H, O, CO system reduces to the O, H sys- 
tem if no CO is allowed to accrete onto grains and the 
OH radical is treated as unreactive. Biham et al. (2001) 
used the notation Pn(i) and Potf) f° r the (independent) 
probabilities that i H atoms and j O atoms are on the 
surface, respectively. For X,Y £ {H, O}, the equation for 
the probability that species X has i atoms is 



^f(*)=fca 



: (X)ri(X)[i\(i-l)-ft(i)] 



+*^ P (X)[(i + l)Px(i + 1) - iPx(i}} 



«x,x 



(±m±H p x (i + 2 ) _ 



(10) 



3.2. Implementation 

The calculations were performed with the use of a Gear 
algorithm on a Cray SV1 computer. To enhance the per- 
formance, the equations were supplied by a subroutine 
which wrote them in an explicit form. Calculations were 
virtually all performed to a time of 1000 yr, which is rather 
short by astronomical scales, but is more than sufficient to 
allow the minor species to reach a steady-state condition 
and the major species to increase linearly with time. We 
have investigated how the needed computer time depends 
on the astronomical time scale and find that for the sys- 
tem studied, there is hardly any difference in computer 
time if the astronomical time is increased 100-fold. To in- 
crease the stability of the calculation, after every call of the 
equation-solving routine, the total sum of the many-body 
probabilities was evaluated and the amount by which it 
deviated from unity was added to the probability for the 
state with no minor species. 



4. Results 

Unless we state that slow H or slow rates are being used, 
the results below are for the fast rates discussed in Sect. 2. 
The fast diffusion rates are emphasized both because they 
allow the rapid production of methanol and because they 
cause a larger discrepancy between the results of the 
standard rate equations and more exact methods. Unless 
stated to the contrary, the temperature is fixed at 10 K. 



+k x .Y(N Y )[(t + l)Px(i + 1) - iPx(t)} 

with initial conditions: Ph(0)=1 and Po(0)=l. Note that 
species X depends on species Y only through its average 
abundance. 

The rate equations for the major species are in the 
form 



d(N, 



OH 



(It 



fco,H x ^ *i p h(h) x hPoih) (11) 
fco.H x (N H )(N ) , 



fcH,H X Y 



h(ii - 1) 



dt ^ 2 

= &H,H X 



ii=2 

(7V H (iVH - 1)) 



= fco.O X 7, P0VI2) 



dt 



(12) 



(13) 



fco,o x 



i 2 =2 

{No(N - 1)) 



For the production of OH, the average abundances of O 
and H appear as independent products rather than as a 
correlation, which would be the case if the minor (H, O) 
species were determined with a two-body probability func- 



4.1. Check of the Independence Approximation 

Before proceeding to our main results, it is interesting to 
check the validity of the suggestion by Biham et al. ( 2001 ) 
that independent probabilities be utilized. We have used 
the simple 0,H system for a comparison among the fol- 
lowing five methods: the Monte Carlo approach (Charnley 
2001, Appendix), the master equation approach, the ap- 
proximation to the master equation approach of Biham 
et al. H2001[) . the rate equation approach, and the mod- 
ified rate equation approach (Stantcheva et al. I2001[) . 
Calculations were performed at a fixed concentration of 
gas-phase atomic oxygen (1 cm -3 ) and a variable concen- 
tration of gas-phase atomic hydrogen. No desorption of 
the products was allowed. In the master equation calcula- 
tions, the maximum allowed numbers {A/i} for H and O 
were never larger than 5. 

Figs.Q] [21 and [3] show the mole fractions of the three 
diatomic molecules, calculated by the various approaches 
as functions of the gas-phase atomic hydrogen abundance. 
For the range of conditions investigated, the exact mas- 
ter equation (Exact ME) and the Monte Carlo simulation 
(Appendix) show excellent agreement for all three species. 
The independent probability approximation of Biham et 
al. H2()()l|) (Approx. ME) follows the inaccurate results of 
the rate equations for O2, is in good agreement with the 
exact approaches for H2 , and is in tolerable agreement for 
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Table 5. Abundances and correlations for the H, O sys- 
tem 



Average 


Exact 


Approx. 


(Nh> 


6.24(-03) 


5.69(-03) 


(No) 


1.82(-01) 


1.23(-08) 


(N H )(N ) 


1.14(-03) 


7.00(-ll) 


<N H No) 


4.80(-ll) 





T 




-2-10 1 2 3 

Log n(H) /cm 3 



Fig. 1. H,0 system. Mole fraction of surface molecular 
oxygen determined via various methods plotted vs the gas- 
phase concentration of H for a 10 K system. 

it approaches the inaccurate rate equation results. In gen- 
eral, the semi-empirical modified rate equation approach 
outperforms the independent probability approximation. 
Biham (private communication) has reported better suc- 
cess with the approximation for slower diffusion rates, so 
one should not rule it out for all situations. 

Another test of whether or not an approximation based 
on independent probabilities is useful is to compare vari- 
ous correlations and averages computed with both the ex- 
act and the approximation master equation approaches. 
In Table [SJ we list the average surface abundances and 
correlations for the atomic H and O when n(H) = 77.(0) = 
1 cm~ 3 . We can see that the values calculated with both 
approaches are in mixed agreement, just as Figs. 1-3 sug- 
gest. For example, while the hydrogen surface abundance 
is almost identical in both cases, the oxygen abundances 
differ from each other by almost seven orders of magni- 
tude. Moreover, the correlation (NhNo) is nowhere near 
the product of the individual averages when calculated by 
the exact approach. 

4.2. H,0,CO system 

For this system, we performed calculations at the three 
densities shown in Table 2] using the simple rate equation, 



£0.5 









JP O-G Monte Carlo 




/ A— A Rate Equations 




ffi y— V Biham et al. 




J □— □ Modified Rates 




■fir x— X Green etal. _ 



-2-10 12 3 4 



Log n(H) /cm 

Fig. 2. H,0 system. Mole fraction of surface molecular 
hydrogen determined via various methods plotted vs the 
gas-phase concentration of H for a 10 K system. 



1 


- ¥-^~=, 0-0 MotlteCarI ° 




/A A— A Rate Equations 




T \ w y'vVu V - V Biham et al. 


c 


/ nr A 13—0 Modified Rates 


_o 




CJ 


vft M\ V Greenetal - 


cS 




£0.5 








o 




% 










-2-10 12 3 4 



Log n(H) /cm 

Fig. 3. H,0 system. Mole fraction of surface OH deter- 
mined via various methods plotted vs the gas-phase con- 
centration of H for a 10 K system. 

ter equation calculations, the set of minimum A/i's which 
must be used is {2,2,1,1,1} , since it is necessary to con- 
sider at least two atoms of H and O on the surface for the 
production of H2 and O2, respectively. Calculations were 
first performed with this minimal cutoff for the five minor 
species, and the results checked by comparison with the 
Monte Carlo method, and by increasing the A/i's. In gen- 
eral, one can get a picture of what cutoff is needed for each 
species by looking at the average number of each minor 
species calculated at a given cutoff, or even the number 
obtained via the rate equation method. If this number ap- 
proaches unity, then a higher cut-off is needed, and if it 
exceeds unity, then it is reasonable to treat the species 
as a major one not requiring inclusion in the many-body 
probability P. 

Because the steady-state criterion leads to the fact 
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density cases are equal, we raised the cut-off of OH from 
1 to 2 leading to A/j's of {2,2,2,1,1}, which should pro- 
duce more accurate results. In the case of high density, 
however, the O and OH abundances are sufficiently high 
(see the detailed discussion below) that we also tried the 
case {2,3,3,1,1}. Any further increase of A/i's led to an 
increase in the computing time without changing the re- 
sults significantly. After the following discussion, we turn 
to analogous calculations with slow diffusion rates, where 
the cutoff problem is more severe. 

The calculated surface populations at 10 K for low, 
medium, and high density are shown, respectively, in 
Tables El El and |H1 for a time of 10 3 yr. In addition to 
the individual populations, the total number of surface 
species is shown, as is the CPU time utilized for the cal- 
culation. It can be seen that at most one monolayer is 
built up during the time of the calculation. With the nor- 
mal assumptions, one monolayer of material corresponds 
to a fractional abundance with respect to the total gas 
density of 10~ 6 . It is to be remembered that only major 
species can be detected on grain surfaces (or in subse- 
quent evaporation into the gas during star formation) so 
that inaccuracies in minor species need not pose a critical 
problem. 

With the large diffusion rates, the minor species (H, 
O, OH, HCO, and CH3O) all have exceedingly low abun- 
dances (< 1 per grain), a situation known as the "ac- 
cretion limit." Such abundances are not capable of being 
determined accurately by the Monte Carlo method, which 
yields integers only. Occasionally we have run the sim- 
ulations several times and averaged the results, so that 
non-integers can be obtained. In the accretion limit, the 
rate equations should not reliably yield accurate answers 
for the stable species, and it is easy to see from the tables 
that such can be the case here. Concerning the more accu- 
rate methods, the minimal master equation method (des- 
ignated 22111; note that the commas have been removed) 
takes about as little CPU time as the rate equation ap- 
proach and never yields results more than a factor of two 
different from the Monte Carlo approach for species which 
can be compared. The latter method is significantly more 
computer intensive. The master equation calculation with 
cutoffs 22211 is typically even better, but at the expense 
of a factor of three in computer time. For the high density 
case, the oxygen atom abundance is 0.5, and increasing the 
cutoffs to 23311 improves the agreement with the Monte 
Carlo method but, again, at the expense of increased com- 
puter time. We conclude that under these physical condi- 
tions, the minimal master equation method is a fast and 
reliable appoach for the H, O, CO system. Although the 
results for the modified rate method were not presented, 
this approach typically does better than the simple rate 
method but worse than the other approaches. 

The tabulated results are all for a temperature of 10 
K. We have also done calculations for high density at 
temperatures through 20 K. The results are plotted for 
major species in Figs. 0] and [S] in terms of mole frac- 




Temperature (K) 

Fig. 4. The mole fractions of surface O2, H2O, and CO 
after 1000 yr for high density conditions plotted vs. tem- 
perature (K) 



Although the major species have surface abundances that 
increase linearly with time, the mole fractions remain con- 
stant. Results are shown for the Monte Carlo (MC), mas- 
ter equation (ME), simple rate equation (RE), and mod- 
ified rate equation (MR) approaches. An analogous plot 
based on the Monte Carlo method is shown in Caselli et 
al. H2002fl and is in a good agreement with the present 
results. The Monte Carlo and master equation results are 
essentially identical, while the rate equation results are at 
best in mixed agreement with the two exact approaches. 
The modified rate method is significantly better than the 
simple rate approach but still can show factor of 2 or 
greater disagreement with the Monte Carlo and master 
equation methods. 

Generally speaking, the results show that under the 
high density conditions considered here, the production of 
methanol is reasonably efficient only at the lower temper- 
atures considered. Its efficiency at 10 K appears to peak 
at intermediate densities (see Table EJ- The mole fraction 
of CO2 also decreases strongly at the higher temperatures, 
while O2 and CO tend generally to increase with increas- 
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Table 6. Calculated populations for surface species at low density and 10 3 yr 



Species 


Rate eq. 


Monte Carlo 


Master Eq. 


Master Eq. 








22111 


22211 


Total 


l.o9(+04) 


l.o5(+04) 


1 GO ( 1 n A \ 

l.o3(+04) 


1.d5(+04) 


lotal (monolayers) 


1.68(-2) 


1.65(-02) 


1.d3(-02) 


l.o5(-02J 


TT 


r.zt^-Uo ) 


t.utn+uu ) 


/.9D^-UoJ 


1 .9D(-Uo ) 


o 


5.21(-07) 


0.00(+00) 


1.92(-02) 


1.90(-02) 


OH 


5.21(-07) 


0.00(+00) 


1.86(-02) 


1.93(-02) 


H 2 


l.ll(+02) 


2.00(+00) 


1.94(+00) 


1.94(+00) 


o 2 


3.64(-07) 


1.70(+02) 


1.62(+02) 


1.62(+02) 


H 2 


1.03(+04) 


9.90(+03) 


9.65(+03) 


9.86(+03) 


CO 


1.51(+02) 


0.00(+00) 


2.81(-02) 


2.81(-02) 


HCO 


3.28(-07) 


0.00(+00) 


1.22(-02) 


1.22(-02) 


H 2 CO 


1.55(+02) 


0.00(+00) 


2.82(-02) 


2.83(-02) 


H 3 CO 


3.28(-07) 


0.00(+00) 


1.23(-02) 


1.23(-02) 


CH3OH 


6.17(+03) 


6.40(+03) 


6.28(+03) 


6.39(+03) 


C0 2 


2.24(-07) 


9.00(+01) 


8.95(+01) 


8.95(+01) 


CPU (s) 


0.5 


11 


1 


3 



Table 7. Calculated populations for surface species at intermediate density and 10 3 yr 



Species 


Rate eq. 


Monte Carlo 


Master Eq. 


Master Eq. 








22111 


22211 


Total 


1.50(+05) 


1.34(+05) 


1.25(+05) 


1.33(+05) 


Total (monolayers) 


1.50(-1) 


1.34(-01) 


1.25(-01) 


1.33(-01) 


H 


5.41 (-06) 


1.00(+00) 


3.01(-03) 


2.88(-03) 


O 


9.75(-06) 


0.00(+00) 


1.35(-01) 


1.36(-01) 


OH 


9.75(-06) 


0.00(+00) 


l.ll(-Ol) 


1.35(-01) 


H 2 


2.20(+01) 


5.00(-01) 


7.32(-01) 


7.01(-01) 


O2 


1.1 2 (-04) 


9.40(+03) 


8.92(+03) 


9.03(+03) 


H 2 


8.57(+04) 


6.02(+04) 


5.18(+04) 


5.93(+04) 


CO 


3.39(+03) 


1.00(+00) 


7.43(-01) 


7.76(-01) 


HCO 


7.36(-06) 


0.00(+00) 


1.14(-01) 


1.14(-01) 


H 2 CO 


3.47(+03) 


1.00(+00) 


6.18(-01) 


7.11 (-01) 


H3CO 


7.36(-06) 


0.00(+00) 


1.21(-01) 


1.22(-01) 


CH3OH 


5.79(+04) 


5.79(+04) 


5.79(+04) 


5.81(+04) 


C0 2 


8.23(-05) 


6.60(+03) 


6.64(+03) 


6.64(+03) 


CPU (s) 


1 


13 


1 


3 



4.3. Slow Diffusion Rates for the H,0,CO system 

As the diffusion rates are lowered, the abundances of sur- 
face oxygen atoms and OH radicals increase dramatically 
and the need for a more detailed treatment than the rate 
equation approach lessens. Indeed, as the average abun- 
dance of a minor species becomes greater than unity, it be- 
comes difficult for the many-body master equation treat- 
ment to converge, so that it is best to remove those species 
from the many-body probability Eq. and compute 
them via equations similar to Eq. © • In Tables El and 
ITUlwe plot results at low and high densities obtained with 
the "slow" (M2) diffusion rates of Ruffle & Herbst 
using the rate equation, Monte Carlo, and master equa- 
tion methods. For low density, we report results for both 
the minimal cutoff case and for a calculation, designated 
2xxll in which O and OH are not treated probabilistically. 



tion since the surface abundances for O and OH are very 
large. 

For the slow rates, it can be seen that little chemistry 
occurs except for the formation of water and, at high den- 
sities, OH. (Any H2 formed has desorbed.) Nevertheless, 
the simple rate equation approach agrees surprisingly well 
with the 2xxll reduced master equation approach. This 
occurs despite the fact that the average H atom abundance 
is less than unity, presumably because the abundances of 
its reactive partners O and OH exceed unity. The Monte 
Carlo method is also in fine agreement for major species 
although it once again is more computer intensive; it can- 
not really be compared with the other methods for minor 
species given the large statistical uncertainties in the small 
numbers. For low density, the minimal cutoff master equa- 
tion approach is unreliable, as is to be expected when the 
O abundance exceeds unity. For high density, the O atom 
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Table 8. Calculated populations for surface species at high density and 10 3 yr 



Species 


Rate eq. 


Monte Carlo 


Master Eq. 


Master Eq. 


Master eq. 








22111 


22211 


23311 


Total 


1.17(+06) 


l.ll(+06) 


9.20(405) 


9.63(405) 


1.09(406) 


Total (monolayers) 


1.17(400) 


1.11(400) 


9.20(-01) 


9.63(-01) 


1.09(400) 


H 


3.42(-10) 


0.00(400) 


1.30(-08) 


1.09(-08) 


8.29(-09) 


u 


4.0Z(-U1 ) 


i.um+uu ) 


a at\( m ^ 
4.911^-111 ) 


£ of; / ni \ 
O.UO(-Ul j 


0. 1 D(-U1 ) 


OH 


4.52(-01) 


1.00(400) 


2.60(-01) 


4.05(-01) 


5.97(-01) 


H 2 


8. 79 (-8) 


0.00(400) 


2.90(-06) 


2.50(-06) 


1.89(-06) 


2 


2.74(405) 


2.81(405) 


1.80(405) 


1.91(405) 


2.68(405) 


H 2 


2.51(405) 


1.79(405) 


1.00(405) 


1.33(405) 


1.71(405) 


CO 


6.47(405) 


5.28(405) 


4.80(405) 


4.93(405) 


5.23(405) 


HCO 


6.71(-04) 


0.00(400) 


2.10(-01) 


1.88(-01) 


1.53(-01) 


H 2 CO 


1.86(402) 


5.01(4-04) 


6.40(404) 


5.98(404) 


5.12(404) 


H3CO 


3.96(-07) 


0.00(4-00) 


5.30(-02) 


4.54(-02) 


3.62(-02) 


CH3OH 


7.33(-02) 


1.10(404) 


1.90(404) 


1.56(404) 


1.17(404) 


C0 2 


2.04(402) 


5.82(404) 


7.90(404) 


7.15(404) 


6.01(404) 


CPU (s) 


1 


30 


1 


3 


10 



Table 9. Calculated populations for surface species with "slow" diffusion rates at low density and 10 3 yr 



Species 


Rate eq. 


Monte Carlo 


Master Eq. 


Master Eq. 








22111 


2xxll 


Total 


1.68(404) 


1.69(404) 


6.98(403) 


1.68(404) 


Total (monolayers) 


1.68(-02) 


1.69(-02) 


6.98(-03) 


1.68(-02) 


H 


8.22(-02) 


0.00(400) 


8.43(-02) 


8.15(-02) 


O 


3.75(400) 


3.00(400) 


7.86(-01) 


3.79(400) 


OH 


3.75(400) 


5.00(400) 


1.81(-01) 


3.79(400) 


H 2 


1.71(-07) 


0.00(400) 


1.60(-07) 


1.50(-07) 


O2 


8.84(-10) 


0.00(400) 


2.66(-ll) 


9.00(-10) 


H 2 


1.03(404) 


1.03(404) 


5.05(402) 


1.03(404) 


CO 


6.47(403) 


6.59(403) 


6.47(403) 


6.47(403) 


HCO 


1.15(-04) 


0.00(400) 


3.93(-05) 


7.38(-05) 


H 2 CO 


1.58(-01) 


0.00(400) 


5.49(-02) 


9.96(-02) 


H3CO 


2.76(-09) 


0.00(400) 


3.27(-10) 


l.ll(-09) 


CH3OH 


2.52(-06) 


0.00(400) 


3.04(-07) 


1.00(-06) 


C0 2 


1.36(-14) 


0.00(400) 


1.24(-15) 


8.77(-15) 


CPU (s) 


0.4 


12 


1.2 


0.2 



significance should include the Eley-Rideal mechanism of 
surface chemistry, in which gas-phase species collide reac- 
tively with nearly stationary species on grain surfaces. 

5. Discussion 

We have shown that a moderately complex network of 
chemical reactions that occur diffusively on the surfaces 
of dust particles can be studied successfully by a master 
equation approach previously used only for significantly 
simpler systems (Biham et al. |2001 Green et al. I2()()l|l . It 
is important to use an "exact" method such as the mas- 
ter equation approach or its Monte Carlo realization when 
there is a very low surface abundance of reactive species, 
since the rate equation method may be inaccurate, and the 
semi-empirical modified rate equation method may not be 
entirely correct either. The advantage of the master equa- 



mer involves the solution of simultaneous differential equa- 
tions. It is therefore facile to consider the gas-phase chem- 
istry occurring simultaneously, because gas-phase abun- 
dances are also determined by solving simultaneous dif- 
ferential equations. Moreover, it is possible to determine 
when the master equation method must be used for all 
minor species, and when it is acceptable to use the sim- 
pler and faster rate equation approach to diffusive surface 
chemistry. Specifically, when the abundances of reactive 
species on grain surfaces begin to exceed unity, there may 
be no need to use the more detailed approach, as can eas- 
ily be tested. Calculations reported here show instances 
where even if the average atomic hydrogen surface abun- 
dance is below unity, the simple rate equation method 
works well if the O and OH abundances are high. Finally, 
for the system studied here, the master equation approach 
is actually faster than its Monte Carlo analog. 
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Table 10. Calculated populations for surface species with "slow" diffusion rates at high density and 10 3 yr 



Species 


Rate eq. 


Monte Carlo 


Master Eq. 
2xxll 


Total 


1.45(+06) 


1.45(+06) 


1.45(+06) 


Total (monolayers) 


1.45(+00) 


1.45(+00) 


1.45(+00) 


H 


2.33(-05) 


0.00(+00) 


2.33(-05) 


O 


4.50(+05) 


4.51(+05) 


4.50(+05) 


OH 


1.97(+05) 


1.97(+05) 


1.97(+05) 


H 2 


1.38(-14) 


0.00(+00) 


1.38(-14) 


2 


4.25(+00) 


4.00(+00) 


4.25(+00) 


H 2 


1.52(+05) 


1.52(+05) 


1.52(+05) 


CO 


6.47(+05) 


6.48(+05) 


6.47(+05) 


HCO 


5.04(-03) 


0.00(+00) 


5.01 (-03) 


H 2 CO 


3.88(-03) 


0.00(+00) 


3.87(-03) 


H3CO 


2.94(-ll) 


0.00(+00) 


2.93(-ll) 


CH3OH 


2.19(-11) 


0.00(+00) 


2.17(-11) 


C0 2 


4.75(-08) 


0.00(+00) 


4.73(-08) 


CPU (s) 


0.6 


27 


0.6 




10 12 14 16 18 20 
Temperature (K) 



Fig. 5. The mole fractions of surface H2CO, CH3OH, and 
CO2 after 1000 yr for high density conditions plotted vs. 
temperature (K) 



With all of these advantages, one might conclude that 
it should be facile to implement the master equation ap- 
proach in current complex gas-grain models of interstellar 
clouds. But this optimism is misplaced. If the master equa- 



still larger systems of reactions, such as that used in cur- 
rent gas-grain chemical models (Ruffle & Herbst 2000), 
some method must be found to reduce the number of si- 
multaneous equations neccessary. Let us consider the ex- 
tent of the problem. The minimal-cutoff (22111) master- 
equation approach for fast diffusion rates requires the si- 
multaneous solution of 3x3x2x2x2 = 72 simultaneous 
differential equations to fully determine the many-body P, 
excluding the additional coupled equations for the major 
species. As one increases the number of minor species, the 
number of simultaneous equations to be solved increases 
dramatically. Suppose, we wish to consider the chemistry 
of deuterium fractionation in the H, O, CO network dis- 
cussed here. Such an extension requires the following new 
minor species: D, OD, DCO, H 2 DCO, HD 2 CO, D 3 CO, 
making a total of 11 such species. The deuterium fraction- 
ation in this model was treated successfully by the Monte 
Carlo method (Caselli et al. 2()()2;. Assuming that mini- 
mal cutoffs of 1 are needed for all these additional species 
with the exception of D, we calculate that a total number 
of 6,912 simultaneous equations is needed for the many- 
body probability function P. This compares unfavorably 
with the total of 652 equations used for both the gas-phase 
and surface chemistry in our most complex models. It is 
difficult to even load the variables for such a calculation 
onto most computers, and the computer time necessary is 
virtually prohibitive. Even if only D is treated stochasti- 
cally and OD, DCO, etc. are treated via rate equations, 
the computer time increases by a factor larger than 10 
compared with the H,0,CO system. 

Given the importance of developing approximation 
methods, we have investigated a simple such approxima- 
tion: the idea that the many-body probability can be 
approximated as the product of individual, independent 
probabilities. We have seen that, even for the simple 0,H 
system, the approximation of independent probabilities 
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body master equation if one assumes fast diffusion rates. 
Although we did not report the results here, our exten- 
sion of the independent probability approach to the more 
complex H,0,CO system has also met with failure. We 
are currently studying other approximation methods. One 
rather promising approach at this time is to limit the total 
number of equations by limiting the total number of reac- 
tive species on a grain surface. We hope to report results 
with this method in the near future. 
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Foundation (US). V. I. S. acknowledges support from grant 
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Appendix A: 

A.l. Kinetics of Chemically Reacting System 

Let us consider an interstellar gas consisting of K dif- 
ferent atomic and molecular species labelled by a^i = 
1 , . . . , K > 1 in a fixed physical volume V with concen- 
tration n ai . Each particle of species is characterized 
by mass m„ velocity Cj, and a set Zj of internal quan- 
tum states. The chemical species in the gas can interact 
through m = 1, M > 1 chemical reactions of the type: 

m : Oi(ci,Zi) +Oij(cj,Zj) -> a fc (c fe ,z fe ) + a ; (c;,z ; ). (A.l) 

The probability that reaction l|A.l|) occurs at a specific 
relative velocity with products scattered in a certain di- 
rection is related to 



dfl 



(A.2) 



where d<r m is the reactive differential scattering cross- 
section, <7y = | c.j — Cj| is a relative velocity, and Q is 
a solid scattering angle. The differential cross section de- 
pends on the interaction potential of the particles involved 
and can be calculated by methods of quantum mechanics 
or measured in laboratory experiments (Light et al. 0-969). 

In astrochemical problems, interstellar gases are usu- 
ally considered under the assumption of local thermal 
equilibrium. This requires that reactive collisions occur 
less frequently than elastic and inelastic processes so that 
a temperature can be maintained. At thermal equilibrium, 
the distribution of possible molecular speeds is given in 
three dimensions by the Maxwellian distribution function 
fL M K It is then possible to obtain the conventional rate 
equations of chemical kinetics (Light et al. [1969): 



[^CKk ^'C*i Q&k , a i ^cti ^OLj Qcti ,atj \ 



(A.3) 



where the functions q., are kinetic rate coefficients for the 
forward (with cross section indicated by — ►) and backward 
(with cross section indicated by <— ) directions of chemical 
reaction lA.il 



Qa k ,a, = J dek dcigki da n 



A.2. Stochastic Approach to a Chemically Reacting 
System 

Chemical kinetics in a rarefied interstellar gas can be 
formulated as a stochastic evolution of an ensemble of 
atoms and molecules. A stochastic approach is based on 
the relationship between two basic ways of describing the 
chemically reacting and evolving gas: (i) the Liouville 
dynamic equation and corresponding kinetic equations 
(Smith 1969) and (ii) the stochastic laws, describing a 
random process and its stochastic Kolmogorov equation 
equivalent (Gillespie ll976l Marov et al. H997J> . For a space- 
uniform gas, the changes in the state of the gas caused 
by instant collisions can be considered to be jump-like 
Markovian processes (Marov et al. 11997(1 , after which the 
state of the system does not contain the memory of how 
the state was reached. 

In a stochastic treatment, the evolution of the react- 
ing system is governed by the so-called chemical master 
equation (Gillespie 119821 van Kampen 1992): 



dt 



M 



= ]T [a m (N™)<XN"\ t)-Om(N)^(N, t)} , (A.5) 



which is linear with respect to the probability density 
</>(N, t) that a system is described by state N at time t. 
The parameters in the equation are defined below: 



a) The state of the system is characterized by 
N(t) = {N 1 (t),...,N K (t)} , 



(A.6) 



where the atomic and molecular populations Ni(t) for 
each species i are random integer variables in the con- 
sidered gas volume V at time t; 

b) The state of the system realized after an instant 
change of molecular populations in accordance with 
the stoichiometric scheme of reaction m is given by 



{■ 



,Ni-l,..., 
...,N k + l,. 



,JV, + 1,...}; 



(A.7) 



c) The probability that a certain reaction m will take 
place in an infinitesimal time interval [t, t + dt] is given 
by the expression a m (JSS)dt, where a m (N) is indepen- 
dent of dt and is equal to: 



a m (N) = V 1 h ij q m . 



(A.8) 



.(Ml , 



(A.4) 



Here, q rn is the rate coefficient for reaction m, and hij is 
a combinatorial factor equal to the number of possible 
pairs of reacting molecules, or NiNj, for 014 ^ ay, and 
N i {Ni-l)/2iora i =a j . 

The function a m (N) depends on the specific chemical 
channel m, the current gas state N = N(t), and the 
gas temperature and volume. This function, usually 
called the propensity function (Gillespie I1976[) . refers 
to processes that lead away from state N. The func- 
tion a m (N m ) depends on the state of the gas N m , and 
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Because N(f) refers to a Mar kovian jump-like random 
process, the time distribution between collisions is given 
by an exponential law: 

F{r(N — > N') < r} = 1 - cxp (-oo(N)t) , (A.9) 
where ao, the total reaction probability, is defined by 

M 

oo(N) = «™( N )- ( A - 10 ) 

m— 1 

Expressions (| A. 6|> - (|A.10|I give an exact definition of 
the random state of the system N(t) describing a chemi- 
cally reacting gas within a stochastic framework. How do 
we describe a system of reactions on the surface of a dust 
particle? 

It is well known that the random nature of interstel- 
lar grain surface chemistry, as well as the accretion and 
desorption processes leading to grain mantle growth, nec- 
essarily involves a stochastic framework (Tielens & Hagen 
115521 Tielens & Charnley 115571 Herbst 125001 Charnley 
2001|. The stochastic treatment of grain surface chem- 
istry can be formulated in terms of the master equation 
approach. To do this, the surface should be represented by 
a lattice or a monolayer - each lattice point corresponds 
to a surface site. A lattice point can assume a number of 
distinct values that stand for various adsorbed molecules 
(with zero for a free site). The monolayer, together with all 
its site populations, is considered as a state vector N for 
the surface reacting system. Each surface reaction changes 
the monolayer population vector in accordance with the 
reaction stoichiometry. The evolution of this surface re- 
acting system over time is described by a chemical master 
equation - Eq. IA.5|I . Since there are a large number of 
reactive sites on interstellar dust particles, additional as- 
sumptions are normally made for stochastic theories - that 
all sites are identical, that one need only follow the num- 
ber of particles of a given species on the entire grain, and 
that the distribution of particles on a grain is random. 

In astrochemical environments at very low tempera- 
tures (~ 10 K), light atoms are the major mobile species 
on grain surfaces. Atoms migrate mainly by thermal hop- 
ping from site to site with a timescale 

r h a =^ 1 exp(Sg/fer d ), 

where E^ is the energy barrier for surface diffusion for 
atomic species a, T<j is the surface temperature, and v 
is the vibrational frequency of the particle in the lattice 
binding site (~ 10 12 s _1 ). Light hydrogen atoms can also 
migrate by quantum mechanical tunneling, with a char- 
acteristic timescale th ~ 10~ 12 s. Migration on the sur- 
face leads to reactions with other light migrating atoms 
or heavy and relatively static molecules and radicals with 
reactive transition probabilities 

a m (N) = x p m x (t- 1 + t^ 1 ), (A.ll) 

where the r _1 factors are for the two reacting species and 
can be either for surface hopping or tunneling. The fac- 
tor p m is unity unless there is an activation energy bar- 



exp(— E m /kTd) or, for the case of quantum mechanical 
tunneling through the potential barrier of height E m and 
width L m , is equal to exp(— 4irL m (2[iE m ) ' 5 /h), where /x 
is the reduced mass. The prime in the combinatorial fac- 
tor h' s „■ means that there is a denominator with the actual 
number of sites on a grain. The denominator converts the 
rate of diffusion from one site to another into the rate of 
diffusion t 7j ff over the equivalent of an entire grain. As sug- 
gested by Charnley |2JJ01|), the surface chemical network 
should be extended by interpreting the accretion and des- 
orption processes as additional reactions responsible for 
the chemical coupling between gas-phase and grain man- 
tle fractions of the interstellar gas. 

A. 3. Methods for Solving the Chemical Master 
Equation 

A. 3.1. Monte Carlo Algorithms 

These algorithms are based on the fact that the probabil- 
ity of generating the stochastic "trajectory" with a Monte 
Carlo algorithm is exactly the probability that would come 
out of the solution of the corresponding master equation. 

To accomplish this, the homogeneous Markovian pro- 
cess N(t) is replaced by an equivalent uniform Markovian 
chain. An exact realization for the Markovian chain on 
a discrete time grid is as follows. We choose a time in- 
terval At and determine times t\ = lAt, <2 = 2At, . . . , 
tg = 6 At, for which we will store the values of Ni = N(ti), 
N 2 = N(i 2 ),- ■ ■ , Na = N(^), respectively. To determine 
the state Na+i of a dust particle, we perform the following 
steps: 

a) we determine randomly: (1) which reaction will 
occur during the transition via (jA~10)) and (jA~TTjl , 
and (2) in accordance with the probability distribu- 
tion (| A. 9|) . the time r that has elapsed from the pre- 
vious transition. This can be done using the direct 
simulation or first reaction methods (Gillespie |l 976 
Lukken et al. lT55%l Charnley 1T55S1 l2"07)T|) . the former 
of which is used here. Then, we update the species' 
populations, and advance the transition time counter 
as TW = T^" 1 ) +t ; 

b) if the transition time counter satisfies the following 
condition: 

ts+i < < t s+2 , 

it means that we have stepped over the next point on 
the time grid. In such a case, we assign a value to Na+i 
that is equal to the current state. If the inequality does 
not hold, we repeat the operations from step a). 

We repeat the whole procedure until the time counter 
reaches a fixed time T /j no / . 

Algorithmic steps (a) - (b) represent the exact Monte- 
Carlo procedure for solving the chemical master equation. 
One realization for the Markovian chain Ni, N 2 , N3, . . . , 
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important that the procedure be repeated and many tra- 
jectory realizations be obtained. Due to the linearity of the 
chemical master Eq. (|A.5ll and, consequently, of analogous 
Monte-Carlo algorithms for its solution, the evolution can 
be calculated by averaging through trajectory realizations 
of the random process N(i). 

A. 3. 2. Direct Solution 

For surface chemistry, one natural way to deal with the 
stochastic approach is to create one probability variable 
for each possible state of the reacting chemical system un- 
der study. If the system can be contained within a limited 
set of possible states so that the number of molecules of 
each species is limited by some fixed value Ni during the 
system evalution - N.^t) < Ni - it is useful to adopt as a 
representation of the state probability distribution 0(N, t) 
a set of many-body probabilities P for specific numbers of 
molecules for the species being considered. By substituting 
this set into the chemical master equation, we can obtain 
a set of coupled differential equations for the time deriva- 
tive of the detailed probabilities of all possible states - 
&P(»i, ...**). 

Such a set of differential equations is used in this pa- 
per, but with a caveat. For surface chemistry, it is nor- 
mally not necessary to include all species in the model in 
the realization of the master equation. Generally, only a 
few reactive species have surface abundances so low that a 
stochastic treament is necessary. All other molecules can 
be described by a deterministic approach with coventional 
rate equations. These two (stochastic and deterministic) 
subsets of the system are coupled through the combinato- 
rial factors 

Ni 

k 3 = (N J )J2*P(-:h-) (A.12) 

i=l 

where (Nj) is a mean molecular population for the de- 
terministic j species and the P(...,i, ...) are probabilities 
for the stochastic I species. This division allows us to com- 
bine the direct solution of the chemical master equation for 
stochastic species with simple or modified rate equations 
for deterministic species to produce realistic grain-surface 
chemical networks. 
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